%the revisible part are K,N,RNK, R(s), sigma(s), y0, tspan
clear all
% parametor of index J
K=10; N=4; RNK=[4,4,ones(1,8)];%,1,1,1,1,1,1,1,1,1,1]; %control index
%K=5; N=4; RNK=[4,4,ones(1,3)];
% K=3; N=4; RNK=[4,2,1];
% construct JNKR
JNKR = index(K,N,RNK);
fname='indexJNKR';
save(fname,'JNKR');
load(fname,'JNKR');
% time 
t = 1;
% construct orthogoal basis in L2[0,t]
orthbasis(t);
% interest
R=@(s) cos(exp(-s));
R=@(s) 1+s-s;
% volatility
alpha=0; m=0; beta=1; sig0=0;

% initial condition
y0=zeros(size(JNKR,1),1);
yic=1;
y0(1)=yic;
% solver
tspan = linspace(0,t,100);
options = odeset('AbsTol',1e-12,'RelTol',1e-6);
%for mean
Y=zeros(100,1);
%for var
Y=zeros(100,size(JNKR,1));

numpath=10;
for ii=1:numpath
    disp(ii)
    sig=ouprocess(alpha,m,sig0,beta,10000);
    [time,y] = ode45(@randomforce,tspan,y0,[],R,sig);
    %Y=Y+y(:,1); %for mean
    Y=Y+y;
end
% error
%meanerr(time',y,R,sig,y0,Y/numpath)
varerr(time',y,R,sig,y0,Y/numpath)


